function [ params ] = calcCellParamsFromNegVstep ( varargin )
% calcCellParamsFromNegVstep
%
% Calculate Rs, Ri, Cm, and tau from a trace with one or more 
% negative-going V steps
%
% % Optional Arguments
%  1: showResponseFig boolean: show/hide the averaged response figure
%  2: currentTrace struct (if not present, will load from gcf.UserData)
%  3: channel number: channel 1 or 2
%
% NB: tailored for specific 4-step pulse with following parameters:
% sampleRate = 10000
% amplitude is extracted; typically -5
% squarePulseTrainNumber extracted, e.g. = 4
% squarePulseTrainISI, e.g. 0.2
% squarePulseTrainDelay, e.g. 0.1
% squarePulseTrainWidth, e.g. 0.1
% Also note that the last 20 msec are assumed to be usable as baseline!

% gs 2008 03 08
% Ben 2009-09-21 display holding current (aka baseline) with cell params
% Ben 2009-09-21 return cell parameters from function
% Ben 2009-09-21 optionally accept currentTrace as argument
% Ben 2009-09-21 optionally hide the averaged response figure
% Taro 2010-02-02 options for channel num
% ----------------------------------

% set this flag to 1 to override faulty stim trace with standard params
% useStandardParams = 1;
useStandardParams = 0;  % BS 2009-07-08 read params from header (my ISI is sometimes 0.25s)

% If no currentTrace argument, load from UserData of current figure
if ( nargin < 1 )
    showResponseFig = true;
    data = get(gcf, 'UserData');
    currentTrace = data.currentTrace;
else
    showResponseFig = varargin{1};
    currentTrace = varargin{2};
end

% trace data
if ( nargin < 3)
    traw = currentTrace.data.ephys.trace_1;
else
    traw = currentTrace.data.ephys.(['trace_',num2str(varargin{3})]);
end
% TODO: add some filtering options

% stimulus parameters
p = currentTrace.header.ephys.ephys.pulseParameters{:};
if isempty(p) % KLUGE!! why is the first item empty??
    p = currentTrace.header.ephys.ephys.pulseParameters{2};
end

sr = p.sampleRate;
V = p.amplitude;
numsteps = p.squarePulseTrainNumber;
isi = p.squarePulseTrainISI;
delay = p.squarePulseTrainDelay;
dur = p.squarePulseTrainWidth;

if useStandardParams
    V = -5;
    numsteps = 4;
    isi = 0.2;
    delay = 0.1;
    dur = 0.1;
end

% calculate baseline and baseline-subtract the whole trace
snippetAtStart = traw(1 : delay * sr - 1);
snippetAtEnd = traw((end - round(0.02 * sr)) : end); % last 20 msec
baseline = mean([snippetAtStart; snippetAtEnd]);
tbs = traw - baseline;

% calculate an average response
for n = 1 : numsteps
    startInd = round((delay + (n-1) * isi) * sr);
    endInd = startInd + round(dur * sr);
    tracelet(:,n) = tbs(startInd : endInd);
end
tmeanbs = mean(tracelet,2);

% make a corresponding time trace
x = 0:length(tmeanbs)-1;
x = x / sr;

% find peak, within first 10 msec
[Ipeak, peakX] = min(tmeanbs(1:round(0.01*sr)));
% assumes pulse is hyperpolarizing; response is downward

% steady state: average over last 20 msec of response
steadystate = mean(tmeanbs(end-round(0.02*sr):end));

% calculate Rseries and Rinput
Rs = 1000 * V / Ipeak;
Ri = 1000 * V / steadystate - Rs;

% steadystate-subtracted mean tracelet:
tmeanss = tmeanbs - steadystate;

% calculate Cm, based on time for signal to fall to 1/e of peak
efold = find(tmeanss(peakX:end) >= (Ipeak-steadystate)/exp(1));
dropInd = efold(1);
tau = x(dropInd);
Cm = 10^6 * (Rs + Ri) * tau / (Rs * Ri);

% graphics
if showResponseFig
    figure('Color', 'w');
    subplot(2,1,1)
    plot(x,tmeanbs, 'b-')
    title('Average current response');
    
    subplot(2,1,2)
    plot(x(peakX:end), tmeanss(peakX:end));
    hold on
    plot(tau, tmeanss(dropInd),'ro');
    title('Average response, baseline steady-state subtracted; 1/e intercept in red');
end

% command line display
disp(' ');
disp('Rs   Ri   Cm   tau   I_hold');
disp([num2str(round(Rs)) '  ' num2str(round(Ri)) '  ' ...
num2str(round(Cm)) '  ' num2str((tau)) '   ' num2str(round(baseline))]);

disp(' ');
disp(['Rs = ' Rs]);

% results -> workspace
cellParameters.Rs = Rs;
cellParameters.Ri = Ri;
cellParameters.Cm = Cm;
cellParameters.tau = tau;
cellParameters.I_hold = baseline;
assignin('base', 'cellParameters', cellParameters); 

params = cellParameters;

